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ABSTRACT 

We present N-body simulations of a new class of self-interacting dark matter models, which 
do not violate any astrophysical constraints due to a non-power-law velocity dependence of 
the transfer cross section which is motivated by a Yukawa-like new gauge boson interaction. 
Specifically, we focus on the formation of a Milky Way-like dark matter halo taken from the 
Aquarius project and re-simulate it for a couple of representative cases in the allowed pa- 
rameter space of this new model. We find that for these cases, the main halo only develops 
a small core (~ 1 kpc) followed by a density profile identical to that of the standard cold 
dark matter scenario outside of that radius. Neither the subhalo mass function nor the radial 
number density of subhaloes are altered in these models but there is a significant change in 
the inner density structure of subhaloes resulting in the formation of a large density core. As 
a consequence, the inner circular velocity profiles of the most massive subhaloes differ sig- 
nificantly from the cold dark matter predictions and we demonstrate that they are compatible 
with the observational data of the brightest Milky Way dSphs in such a velocity-dependent 
self-interacting dark matter scenario. Specifically, and contrary to the cold dark matter case, 
there are no subhaloes that are more concentrated than what is inferred from the kinematics 
of the Milky Way dSphs. We conclude that these models offer an interesting alternative to the 
cold dark matter model that can reduce the recently reported tension between the brightest 
Milky Way satellites and the dense subhaloes found in cold dark matter simulations. 
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1 INTRODUCTION 



The collisionless nature and low intrinsic thermal velocities of 
dark matter (DM) particles are two of the fundamental hypothe- 
ses underlying the remarkably successful Cold Dark Matter (CDM) 
paradigm. The latter of these hypotheses is related to a filtering 
mass scale that sets the minimum mass for self-bound DM struc- 
tures (haloes) to be typically of 0(1M®) and hence the CDM 
model predicts the existence of a large population of Earth-mass 
haloes in the local Universe. Although such prediction has yet to 
be unambiguously tested by observations, there is evidence of a po- 
tential tension with the actual number of observed dwarf galaxies, 
both within the Milky-Way (MW) (the so called "missing satellite 
problem", e.g.lKlypin et al .|1999| |Moore et al.|1 999), and in the 
field ([Zavala et al. 2009 ; Papastergis et al.|2011| as inferred from 
the HI velocity function measured by the ALFALFA survey). This 
challenge to the CDM model is undoubtedly connected to the de- 
velopment of a successful theory of galaxy formation at the scale 
of dwarfs. For instance, it has been shown that the combination of 
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astrophysical processes that inhibit star formation are able to solve 
the "missing satellite" problem (e.g. Koposov et al. 2009 1, but the 
tension with dwarf galaxies in the field remains despite attempts to 
explain it without modifying the CDM model ( |Ferrero et al.|201 1) . 
It remains to be seen if galaxy formation models based on N-body 
simulations that seem to provide a good fit to the luminosity func- 
tion at low masses \Guo et al.|20TT) are also successful in reproduc- 
ing the HI observations of the velocity function. Such a prevailing 
challenge strengthens the case for a solution based on a suppres- 
sion of primordial small-scale density fluctuations. For example, 
models where DM is "warm", i.e., where DM particles have larger 
primordial thermal velocities, are able to alleviate the overabun- 
dance problem (e.g. |Bode et al.|2001||Z~avala et al.|2009f and at the 
same time, if their masses are > 1 keV, avoid current constraints 
on the matter power spectrum as measured by the Lyman— a forest 
( |Boyarsky et al.|2009| >. 

The second hypothesis, that of DM being essentially collision- 
less, has also been contested. The idea of self-interacting Dark Mat- 
ter (SIDM) was first suggested by|Carlson et al.| < |T99"2"| >; |Machacek| 
|et al.| (1993); |de Laix et al.| (|1995|> near ly 20 years ago. More 
than a decade ago, Spergel & Steinhardt (2000) realised that the 
aforementioned overabundance problem of dwarf galaxies, and the 
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observation of low surface brightness (LSB) galaxies having den- 
sity cores (de Blok & McGaugh 1997 1 (which contradicts the high 
density cusps predicted by the CDM model ( |Navarro et al.|1996] 
|1997| >) could be avoided if DM would be self-interacting. Current 
observations of LSB galaxies (Kuzio de Naray & Spekkens 201 1 1 
and MW dwarf spheroidals (dSphs) ( Walker & Penarrubia 2011) 
seem to confirm the presence of density cores in low-mass haloes. 
Since these galaxies are DM-dominated, it is challenging to invoke 
baryonic processes as the main mechanisms responsible of alter- 
ing so drastically the inner density profile of haloes. If DM is not 
cold, then haloes are expected to have cores, although those asso- 
ciated with currently allowed WDM models seem to be too small 
to explain the o bserved cores of LSB galaxies (|Kuzio de Naray &| 
|Spekkens|2011||Villaescusa-Navarro & Dalal|201 1) . Haloes within 
WDM simulations, although less concentrated, seem to have simi- 
lar profiles as their CDM counterparts and they do not show a clear 
sign of developing cores (e.g. |Coh'n e t al. 2000). High resolution 
simulations within the WDM cosmogony are however challenging 
( |Wang & White|20 07) and hence a consensus on the inner density 
profile of WDM haloes has not been reached. SIDM models on the 
other hand, lead naturally to the development of a substantial core 
as was already shown by the first SIDM simulations ( jYoshida et al.| 
|2000a|b||Dave et al.|2001 [|Colm et al.|2002> . It remains to be seen if 
SIDM models are able to explain the observed cores of MW dSphs 
and LSB galaxies. 

The first SIDM models assumed a constant scattering cross 
section and were quickly abandoned since those that could solve the 
small-scale CDM problems seemed to violate several astrophysical 
constraints, such as the observed ellipticity of the mass distribu- 
tion of galaxy clusters (e.g. Miralda-Escude 2002 1 and the surviv- 
ability of satellite haloes (e.g. Gnedin & Ostriker 2001 ). To avoid 
such constraints, simple ad hoc velocity-dependent cross sections 
of the form l/v a were explored (e.g. Col in et al.|2002) , yielding 
encouraging results that however lacked a proper underlying par- 
ticle physics model. It has also been claimed that these velocity- 
dependent SIDM models are not able to solve simultaneously the 
core problem in DM-dominated systems and the "missing satellite 
problem" (e.g. D'Onghia & Burkert 2003 1. 

|Loeb & Weinerj {201 1[ > proposed that the possible existence 
of a Yukawa potential among DM particles can resolve the chal- 
lenges facing SIDM with a constant cross section. The velocity de- 
pendence of scattering through the massive mediator of this dark 
force (similar to a screened Coulomb scattering in a plasma) could 
make scattering important at the low velocity dispersion of dwarf 
galaxies but unimportant at the much higher velocities encountered 
in galaxy clusters. The possibility of exothermic reactions could 
in addition introduce a special velocity scale around which the in- 
fluence of the DM interaction peaks. The existence of dark forces 
was studied earlier as a solution to cosmic ray anomalies through 
enhanced dark matter annihilation ( Arkani-Hamed et al. 2009)). 

A recent analysis by Boylan-Kol chin et al.| ( |20lTa| puts for- 
ward an additional challenge to the CDM model. The authors used 
simulated MW-like haloes in a CDM cosmology to show that the 
observed MW dSphs are only consistent with inhabiting relatively 
low-mass CDM subhaloes, leaving a population of more massive 
subhaloes with no galaxies associated to them, i.e., massive sub- 
haloes of CDM MW-like haloes seem to be too dense to host the 
bright MW dSphs. In a recent extension to this analysis, |Boylan-| 
|Kolchin et al.| ( |2011b[ ) showed that this problem is unlikely to 
be solved invoking standard galaxy formation processes based on 
CDM. This is one of the most serious challenges faced by the CDM 
model and can perhaps be solved by invoking WDM ( Lovell et al.| 



|201 1\ or, alternatively, also naturally avoided in certain SIDM mod- 
els as we explore in this work. 

This is the first paper in a series whose objective is to study the 
properties of DM (sub)haloes within allowed velocity-dependent 
SIDM (vdSIDM) models using state-of-the-art numerical simula- 
tions. In this work we only focus on the recently discovered prob- 
lem pointed out by |Boylan-Kolchin et al.| |201 la|b^ and demon- 
strate how vdSIDM models can reduce the discrepancy between 
observation and theoretical prediction. In Section [2] we briefly de- 
scribe the particle physics model we use and present our numer- 
ical algorithm including test simulations. In Section [3] we present 
simulations of a MW-like dark matter halo using the initial con- 
ditions of one of the Aquarius haloes (Springel et al. 2008} and 
compare its resulting subhalo population, for a couple of represen- 
tative cases in the parameter space of the vdSIDM model, with the 
ones of the standard collisionless CDM model. These models are 
then contrasted against observational data from the MW dSphs in 
Section [4] Finally a summary and the conclusions of our work are 
given in Section[5] 



2 METHODOLOGY 

2.1 Velocity-dependent SIDM (vdSIDM) models 

We use a simplified particle physics model where the self-scattering 
between DM particles of mass m x is set by an attractive Yukawa 
potential with coupling strength a c mediated by a new gauge bo- 
son of mass (either scalar or vector) in the dark sector. We re- 
fer the rea der tojFeng et aT^2010b) ; |Finkbeiner et al!] ( [201 l| l; [Loeb| 
|& Weiner] ( |201 1 1 for details on the particle physics model and its 
available parameter space. If we only consider elastic interactions, 
the scattering problem is analogous to the screened Coulomb scat- 
tering in a plasma, which is well fitted by a transfer cross section 
given by: 



t 7 P 2 ln(l + /?' 1 ) 



P 2 (1 + 1.50 1 



/3 < 0.1 



O.K P < 10 s (1) 



(ln/3 + 1- llrC 1 ?) 2 , /3>10 3 
Lx/« 2 = 2a c m <t> /{m x v 2 ) 



where (3 = 7 

22.7/m0, and v is the relative velocity of the DM particles. Here 
Umax is the velocity at which (cr-ru) peaks at a transfer cross section 
equal to cr™ ax . 

The value of <r™ ax /m x is constrained by different astrophys- 
ical measurements, the most stringent being: i) the observed ellip- 
soidal shape of haloes as implied by X-ray data of galaxy clus- 
ters and ellipticals jMiralda-Escude||2002| |Feng et al.||2010a| >; ii) 
avoidance of the gravofhermal catastrophe leading to inner den- 
sity profiles even steeper than those predicted in CDM ( |Firmani| 
et al. 2001); and iii) avoidance of the destruction of subhaloes 
through collisions with high velocity particles from a larger par- 
ent halo ( Gnedin & Ostriker 200 1}. There is a sum mary of these 
and other constrains in Table I of |Buckley & Fox| ( |2010) and in 
Figure 2 of |Loeb & Weiner (201 1): on the scales of dwarf galaxies, 
&vei ~ 10 km s _1 , the allowed values for the transfer cross section 
are roughly constrained from above by er™ ax /m x < 35 cm 2 g _1 , 
and are much lower at <r ve i ~ 100 km s -1 , where the constraints 
are stronger by approximately two orders of magnitude. Since we 
are interested in the possibility of producing cored density profiles 
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Figure 1. Left panel: Analytic scatter rate profiles for Hernquist haloes of different masses following the average cosmological concentration-mass relation. 
Dashed lines show the result for a velocity-dependent cross section (ti ma x = 30 kms - 1 , <r™ ax /m x = 10 cm 2 g" 1 ), whereas solid lines show the rates 
for a constant cross section (<r™ ax /m x = 10 cm 2 g —1 ). The remaining two panels show a comparison of these scatter rates with those obtained from a 
frozen N-body representation. The simulations were run for 1 Gyr, and the different halo profiles are all offset by one dex for clarity. Solid lines show the 
N-body result, whereas black dashed lines represent the analytic result repeated from the leftmost panel. Middle panel: constant cross section. Right panel: 
velocity-dependent cross section. 



for the haloes associated with the MW dSphs, we will take two 
benchmark points in the (cr™ ax /m x , u max ) parameter space close 
to the aforementioned constraints that maximise the self-interaction 
at the typical velocity dispersion of these dwarfs (see Table [T]in 
Section [3T| below). 

In this work we only consider elastic scattering leaving the 
cases of excited states and their associated exo- and endothermic 
interactions for a future analysis. 

2.2 Numerical technique 

To account for DM self-interactions we follow a standard Monte 
Carlo approach similar to previous implementations ( Kochanek & 
White|2000HBurkert|2000||Yoshida et al.|2000a|b||Dave et al.|2001] 
Craig & Davis|2001||CoIin et al.|2002||D'Onghia et al.|2003[[Koda| 
& Sha piro|2011} , but different from fluid smoothed particle hydro- 
dynamics approaches as in Moore et al. (20001 and Yoshid aet al.| 
P000a| >. 

We determine the scattering probability for every particle i 
with each of its k — 38 ± 5 nearest neighbour^ j in a time step 
AU by 

Pij = — Win-j, hi) a T (vij)vij Ati, (2) 
m x 

where m, is the simulation particle mass, Vij is the relative velocity 
between particles i and j, OT/ra x is the scattering cross section 
per unit mass described in Section [2~T| hi the smoothing length 
enclosing the k nearest neighbours of particle i, and W(rij, hi) — 
w(rij/hi) is the cubic spline Kernel function in 3D normalisation: 

C l-6q 2 +6q 3 , 0<<K |, 
w(q) = -\ 2(l-q) 3 , 5<?<1, (3) 

n { 0, q>l. 

The time step At; is chosen small enough to avoid mul- 
tiple scatterings during a time step by requiring that Ati < 

k (pi (JT {(J vcl 

,i)/m x <T V ci,i) \ where cr V ei,! is the local velocity 

1 This choice is to speed up the neighbour search, but we checked that it 
does not affect any results. 



dispersion at the position of particle i calculated based on its k 
neighbours, and we set k = 10~ 2 , which is sufficiently small to 
avoid multiple scatterings during a step and usually smaller than 
the time step inferred from the dynamical time scale. The total 
probability of a particle to interact with any of its neighbours is 
given by Pi = . Pij /2, where the subscript j is limited to the k 
neighbours. The factor 1/2 accounts for the fact that a scatter event 
always involves two particles, and we therefore need to divide by 
two to reproduce the correct scatter rate. We say that a collision 
takes place between particle i and one of its k nearest neighbours 
j if a; ^ Pi, where a; is a uniformly distributed random number in 
the interval (0, 1). To select the neighbour j that is chosen for colli- 
sion we sort them according to their distance to particle i and select 
the first neighbour I that satisfies x ^ Pij. In the following we 
assume that the self-interaction is isotropic. 

In the case of elastic scattering once a pair is tagged for colli- 
sion we assign to each particle a new velocity given by: 

Vi = V cm + («ij/2) e, 

Vj = Vcm - {vij/2) e, (4) 

where v cm is the centre-of-mass velocity of the pair and e is a unit 
vector that we randomly draw from the unit sphere. This procedure 
conserves energy and linear momentum, but not angular momen- 
tum. We have implemented this numerical scheme in GADGET-3 
(last described in Springel 2005}. 

To test our implementation we apply it first to isolated haloes. 
For a region of volume V, the total number of scattering events is 
given by: 

Ttot = / (otv) (x) dV, (5) 

Jv ^ m x 

where p(x) is the local DM density and (cttv) (x) is the local ther- 
mal average of the transfer cross section times the relative veloc- 
ity. In the non-relativistic limit this is given by an average over a 
Maxwell-Boltzmann distribution function: 

(<j T v) (x) = 3 1 / (cr T v)v 2 e~ v /4tT ™iM d w , (6) 

where cr vc i(x) is the local velocity dispersion. For given density 
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Figure 2. Density (left panels) and velocity dispersion profiles (right panels) of haloes of different masses. The top panels are for the case of a constant 
cross section (<7™ ax /m x = 10 cm 2 g — 1 ) showing the profiles after 25 to. Bottom panels are for the case of a velocity-dependent cross section (i> ma x = 
30 kms -1 , cr™ ax /rri x = 10 cm 2 g — 1 ) after 1 Gyr. In scaled units, the constant cross section curves for all masses collapse to a single one. For the 
velocity-dependent case, evolution progresses faster for lower mass systems, because (<r T ti) peaks at a velocity of 30 km/s. 



and velocity distribution functions we can now calculate the num- 
ber of expected scattering events and compare this to the N-body 
/ Monte Carlo results obtained with the technique presented in the 
paragraphs above. 

As an example of the number of scattering events expected in 
a DM halo, we take a smooth spherical distribution of DM with a 
Hemquist density profile (Hernquist 1990): 

, , Ma 1 

where M is the total mass of the halo and a its scale length. The 
velocity dispersion profile for the Hernquist halo follows from the 
Jeans equation, which for an isotropic velocity distribution and us- 



ing Eq. (J7]l gives: 

GM 
~12a~ 



i(r) = 



Ylrir + a) 



■In 



r + a 



(8) 



r + a 



'25 + 52 (-) +42 (r 



+ 12 



It is then straightforward to compute the scattering rate using 
Eq. |5j. To compare these analytical expectations with N-body 
simulations, it is necessary to take into account the mass resolu- 
tion of the simulation. We therefore need to multiply Eq. l|5} with 
m x /iridm, where mdm is the DM particle mass of the simulation, 
which yields the number of scatter events in the simulation volume. 
The left panel of Figure [T] shows the analytically calcu- 
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lated scatter rate for a constant (<r™ ax /m x = 10 cm 2 g _1 , solid 
lines) and velocity-dependent (w ma x = 30 kms~ , u™ ax /m x 
1 1 ! cm 2 g" 1 , dashed lines) cross section. In the constant cross sec- 
tion case we neglect the velocity dependence, whereas the velocity- 
dependent case assumes the parametrisation of Eq. [T] We show 
these rates for 6 different haloes, where we have chosen the scale 
radii a to follow a cosmol ogically motivat ed mass-concentration 
relation c oc M~ 1,/9 (e.g. Neto et al.|2007 i. In the other two pan- 
els of Figure [T] we compare the analytic scatter rates to the result 
of the N-body (TV = 10 6 , Plummer equivalent softening length 
e = 0.0063 x a) calculation assuming a frozen particle configura- 
tion. The middle panel shows the constant cross section, whereas 
the right panel shows the result for the velocity-dependent cross 
section. All cases show good agreement between the analytic result 
and the simulation. 

In Figure|2](two top panels) we show the density and velocity 
dispersion profiles for the same haloes and the same constant cross 
section with a live N-body (N = 10 B , e = 0.0063 x a) system 
after 25 to, where 



Po 



M 
27ra 3 ' 



vo 



i\/AnGp , t 1 



CTT . / GM 3 

-a\l — =— , (9) 



and a = 2.26 (see |Koda & Shapiro 20111. As expected all haloes 
evolve in the same way if time is expressed in units of to, density in 
units of po and velocity in units of vo- Core expansion is maximised 
around 25 to, and after that the core-collapse phase slowly starts. 
The two lower panels of Figure|2]show density and velocity disper- 
sion profiles for the velocity-dependent case after 1 Gyr. Within 
this time span the haloes reach different levels of core expansion 
due to their different velocity structure with the lowest mass halo 
showing the largest core at that time. 



3 COSMOLOGICAL SIMULATIONS OF A GALACTIC 
SIDM HALO 

3.1 Initial Conditions and Models 

The initial conditions for our cosmological MW-like halo simula- 
tions are taken from the Aquarius Project (Springel et al.||2008] >. 
This has the advantage that we can compare subhalo properties 
object-by-object with CDM simulations that have already been car- 
ried out at very high resolution. In addition, we can directly com- 
pare our results with those of Boylan-Kol chin et al.| ( |201 Tb^ , who 
found the massive subhalo failure problem described in the Intro- 
duction using Aquarius dat^] The Aquarius initial conditions use 
the following cosmological parameters: Q m = 0.25, Qa = 0.75, 
h = 0.73, as = 0.9 and n s — 1; where Q m and Qa are the contri- 
bution from matter and cosmological constant to the mass/energy 
density of the Universe, respectively, h is the dimensionless Hubble 
constant parameter at redshift zero, n s is the spectral index of the 
primordial power spectrum, and as is the rms amplitude of linear 
mass fluctuations in 8 h^ 1 Mpc spheres at redshift zero. The haloes 
were selected to be representative of the MW halo with different 
merger histories. Here we focus on the Aq-A halo and re-simulate 
it within the SIDM model described before at different resolutions. 
We note that we only consider scatter events between high reso- 
lution particles, i.e. massive boundary particles behave like CDM. 



We note that Boylan-Kolchin et al. 



data from the Via-Lactea simulations Diemand et al. 2008 1 and found the 
same problem. 



Name 


Type 


(T™ ax /m x [cm 2 g 1 ] 


%ai [kms x ] 


RefPO 


CDM 


/ 


/ 


RefPl 


SIDM (ruled out) 


10 


/ 


RefP2 


vdSIDM (allowed) 


3.5 


30 


RefP3 


vdSIDM (allowed) 


35 


10 



Table 1. Reference points and their particle physics parameters explored 
in our simulations. RefPl serves only as a benchmark point for tests, since 
it is well-known that such a large constant cross section violates various 
astrophysical constraints. RefP2 and RefP3 do not violate any constraints 
and potentially have a significant effect on the density profiles of low-mass 
subhaloes. The latter two reference points are therefore the ones we will 
mainly focus on. 



This does not introduce any bias effects since our SIDM models 
produce the same large scale structures as CDM, i.e. tidal effects 
due to the boundary particles are treated correctly even when ne- 
glecting their self-scattering. Furthermore the high resolution re- 
gions of all the simulations presented here are free of any contami- 
nation from boundary particles. 

In Table [T] we show the SIDM models we consider in this pa- 
per. RefPO is the vanilla CDM case and RefPl only serves as a test 
case to demonstrate various effects since is ruled out based on astro- 
physical constraints. For instance, observations of the Bullet cluster 
place an upper limit to the DM scattering cross section which is be- 
low the value used in the RefPl case: aT/m x < 1.25 cm 2 g _1 
I Randall et al.|200 8p1 RefP2 and RefP3 on the other hand, do not 
violate any constraints and potentially have a significant effect on 
the density profiles of low-mass subhaloes. The latter two reference 
points are therefore the ones we will mainly focus on. 

In the following, we investigate the effects of DM self- 
interaction on the present-day properties of the MW-like halo and 
its subhaloes, which consist of self-bound groups of particles iden- 
tified by the SUBFIND algorithm ( Springe Tet al.|2001J down to a 
group with a minimum of 20 particles. 



|201 lb) also did the analysis using 



3.2 Main halo 

A first visual impression of the structure of the halo in the differ- 
ent DM models is given in Figure [3] where we show density pro- 
jections of the halo for the models given in Table [T] Clearly, the 
disfavoured RefPl model with a large constant cross section has a 
very different density distribution with a spherical core in the cen- 
tre. Substructures also share these properties being less dense and 
more spherical in this simulation. The vdSIDM models RefP2 and 
RefP3 on the other hand can hardly be distinguished from the CDM 
case (RefPO) in this figure. As we see below, the internal struc- 
ture of the massive subhaloes does actually change significantly 
with respect to the CDM case. We note that the spherical cores in 
self-interacting models are due to the assumed isotropic scattering 
process, which tends to isotropise particle orbits leading to a more 
isotropic and spherical configuration. 

The left panel of Figure|4]shows the main halo density profiles 

3 A large constant cross section model will also create a core and isotropise 
the central part of galaxy clusters. There are stringent constraints on the 
cross section for these systems, for instance, Miralda-Escude 1 2002 1 derived 
a limit of cry / m x < 0.02 cm 2 g — 1 by analysing the ellipticity of a cluster 
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RefPO (CDM) 




(RefPl (S IDM-ruled c 




RefP2 (vdSIDM-allowed) 




RefP3 (vdSIDM-allowed) 




Figure 3. Density projections of the Aq-A halo for the different DM models of Table[T](RefP0-3). The projection cube has a side length of 270 kpc. Clearly, 
the disfavoured RefPl model with a large constant cross section produces a very different density distribution with a spherical core in the centre, contrary to 
the elliptical and cuspy CDM halo. Also, substructures are less dense and more spherical in this simulation. The vdSIDM models RefP2 and RefP3 on the 
other hand can hardly be distinguished from the CDM case (RefPO). 



for the different models, whereas the right panel shows the mean 
free path A = (p (<tt / m-x ) ) ~ 1 as a function of radius for the SIDM 
models. The dotted, dashed and solid lines show different levels 
of resolution, characterised by a particle mass m p and a Plummer 
equivalent gravitational softening length e: Aq-A-5 (m p = 3.143 x 
10 6 M , e = 684.9 pc), Aq-A-4 (m p = 3.929 x 10 5 M , e = 
342.5 pc) and Aq-A-3 (m p = 4.911 x 10 4 M , e = 120.5 pc). The 
runs show good convergence for radii larger than 2.8e indicated by 
the vertical lines. 

In the figure we see that RefPl develops a large core reach- 
ing the solar circle (~ 7 kpc). This is because the cross section 
has no velocity dependence in this case and the particle scattering 



works at full strength irrespective of (sub)halo mass. Although this 
case is ruled out by current astrophysical constraints (see Section 
|2.1[ >, it serves as a reference for the effect of a large scattering cross 
section at the scales of MW-like haloes in a full cosmological sim- 
ulation. On the contrary, RefP2 and RefP3 result in a main halo 
whose density profile follows very closely the one from the CDM 
prediction of RefPO down to 1 kpc from the centre. At smaller radii, 
where the typical particle velocities are smaller, self-interaction is 
large enough to produce a core. The mean free path radial profile 
clearly illustrates the radius where collisions are more important 
for the different SIDM models, which is around the core radius. It 
also highlights the difference between the RefP2 and RefP3 mod- 
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Figure 4. Left panel: Density profile of the Aquarius Aq-A main halo at resolution levels 5 (dotted), 4 (dashed) and 3 (solid) for the different SIDM reference 
points we consider (see Table ^ as shown in the legend. Right panel: Mean free path as a function of radius for the SIDM models we used. The softening 
length (2.8 times Plummer equivalent) of each resolution is marked by a vertical black line. The standard CDM profiles (RefPO) are shown in black. We 
achieve good convergence in all our runs, with the inner profiles changing significantly depending on the DM model employed. Clearly, RefPl produces the 
largest difference having a large core due to the constant scattering cross section. We note that this model is ruled out by current astrophysical constraints and 
is shown here just as reference. 



main halo 




r [kpc] 

Figure 5. Radial profiles of the mean number of scatter events for the Aq- 
A-3 simulations RefPl, RefP2 and RefP3. The large and constant cross- 
section of RefPl produces a significantly higher number of scatter events 
at a given radius compared to RefP2 and RefP3. We note that the vdSIDM 
points RefP2 and RefP3 lead to slightly different scatter profiles. 



els, with the former having a larger core than the latter, because 
its self-interaction cross section peaks at a larger velocity disper- 



sion (occurring at larger radii) despite of having a lower value of 

<jt /rn x . 

Each simulation particle records the total number of scatter 
events during its dynamical evolution. We can use this information 
to show alternatively the size of the collisional radius by construct- 
ing a radial profile of the mean number of scatter events as a func- 
tion of radius. This is shown in Figure|5]for the highest resolution 
Aq-A-3 simulations for RefPl, RefP2 and RefP3. Clearly, the large 
and constant cross-section of RefPl produces a significantly higher 
number of scatter events at a given radius compared to RefP2 and 
RefP3. 



3.3 Subhaloes 

Our main focus in this work is the structural change of the sub- 
halo population in a SIDM halo. In the following we will mainly 
focus on the subhalo population within 300 kpc halocentric dis- 
tance. In the left panel of Figure |6] we first show the ratio of the 
subhalo mass function of our different models (for the highest res- 
olution run, level 3) to the best fit subhalo mass function of the 
CDM Aquarius haloes (see Springel et al. 2008 ). The disfavoured 
RefPl model differs significantly from the CDM prediction, con- 
trary to the RefP2 and RefP3 models that have mass functions very 
similar to the RefPO model. This means that subhalo masses do 
not change significantly in the vdSIDM models, whereas subhaloes 
lose mass in the RefPl model. The right panel of Figure [6] shows 
the subhalo radial number density profiles for the different mod- 
els. The subhalo population has been divided in three mass bins: 
(10 6 , 10 7 ) M Q , (10 7 , 3.2 x 10 8 ) M Q and (3.2 x 10 8 , 10 11 ) M , 
from bottom to top, respectively (the last two bins are shifted up 
by one (two) dex for clarity, in fact the three curves lie on top of 
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Figure 6. Left panel: Ratio of the subhalo mass function for our different SIDM reference models to the best fit subhalo mass function of Spring eTet al.| 
1 2008 1. Right panel: radial distribution of the subhalo number density (the virial radius is marked with a solid black line). The subhalo population has been 
divided in three mass bins: (10 6 , 10 7 ) Mq, (10 7 , 3.2 X 10 8 ) M Q and (3.2 X 10 8 , 10 11 ) M©, from bottom to top, respectively (the last two bins are shifted 
up by one (two) dex for clarity). The constant cross section SIDM model (RefPl) gives a significantly different subhalo abundance, particularly in the inner 
regions of the main halo, than the CDM model (RefPO). On the contrary, the velocity-dependent models (RefP2 and RefP3) are nearly indistinguishable from 
the collisionless case. For clarity, we only include the results for the highest resolution simulations (level 3). 



each other) . Once again, the velocity-dependent scenarios are es- 
sentially indistinguishable from the CDM model (at all masses), 
whereas the RefPl model shows a considerably lower abundance 
of subhaloes within the virial radius of the main halo (particularly 
for the least massive subhaloes), which is more acute in the inner 
regions. Although not shown here, we checked that the lower res- 
olution runs converge to the level 3 curves for both subhalo mass 
function and radial distribution for all models. 

The subhalo evaporation seen in the RefPl model is explained 
by the interaction that takes place between the particles of a given 
subhalo and the higher velocity particles of the surrounding host 
halo. Since these collisions are elastic, the net effect is a transfer 
of energy to the particles in the subhalo, unbinding the ones with 
the lowest binding energies. These interactions are more common 
in the inner region of the main halo where the density is higher, 
and strongly reduce the abundance of subhaloes there. On the other 
hand, evaporation is not effective in the velocity-dependent models 
where self-interactions are strongly suppressed for high relative ve- 
locities, typically expected between particles in the subhaloes and 
those of the host. This is an important result since it implies that 
the overabundance of dwarf haloes is still present for the vdSIDM 
models explored here. This means that the "missing satellite prob- 
lem" would need to be solved by invoking astrophysical processes 
as is done for the CDM model. We note however that this is only 
true for elastic scattering, for the inelastic case, subhalo evapora- 
tion should have a relevant role. We leave the study of this case for 
a future analysis. 

The internal structure of the subhalo population can be appre- 
ciated in the left panel of Figure [7] that shows the density profiles 
at z — for the 15 most massive subhaloes in our highest resolu- 
tion simulations (level 3) for the cases RefPO (black), RefPl (red), 
RefP2 (blue) and RefP3 (green). The figure shows that the vdSIDM 



cases produce subhaloes with cores of approximately 600 pc; there 
is little difference between cases RefP2 and RefP3, with the latter 
having slightly larger cores. Beyond 1 kpc, these cases are indistin- 
guishable from CDM. The case where self-interaction is indepen- 
dent of velocity (RefPl) results in cores that have lower density and 
are more extended (approximately a factor of ~2). The right panel 
of Figure|7]shows the mean free path as a function of radius for the 
difference SIDM cases. We note that one of the shown subhalo of 
the RefPl halo entered already the core-collapse regime, and there- 
fore shows a very steep density profile. Other subhaloes are still in 
the core-expansion phase resulting in significantly shallower pro- 
files compared to the CDM (RefPO) curves. 



We note that the cores of the main halo and the massive sub- 
haloes (left panel of Figure|4]and Figure|7J are roughly of the same 
size. Contrary to previous velocity-dependent SIDM models, ours 
seem to predict no strong scaling of the core size with mass. A con- 
cern might be that for cluster-sized systems, the hierarchical nature 
of structure formation implies that they were assembled from small 
mass halos that formed at earlier epochs. This means it is possible 
that although our vdSIDM models predict a vanishing cross-section 
for cluster scales, such systems might still exhibit large cores due 
to the assembly of smaller significantly cored systems; such a large 
cluster core might violate observational constraints. However,the 
mass contribution from these smaller systems to the central region 
of the cluster is expected to be subdominant ( [Gao & White|2006| > 
and therefore the evolution of the inner regions of larger haloes is 
likely to be affected only in a minor way by the merging of smaller 
cored subhaloes. 
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Figure 7. Subhalo density (left panel) and mean free path (right panel) at z = for the top 15 most massive subhaloes (largest peak circular velocity) 
in our highest resolution simulations (level 3). The vdSIDM models produce subhaloes with cores of approximately 600 pc and there is little difference 
between RefP2 and RefP3, with the latter having slightly larger cores. Beyond 1 kpc, these cases are indistinguishable from CDM profiles. The case where 
self-interaction is independent of velocity (RefPl) results in cores that have lower density and are more extended (approximately a factor of ~2). We note that 
one of the shown subhaloes of RefPl entered already the core-collapse regime resulting in a very steep density profile, whereas most of the other subhaloes 
are all still in the core-expansion phase. 



4 SUBHALO POPULATION: COMPARISON WITH THE 
BRIGHT MW DWARF SPHEROIDALS 

To check the the consistency between the subhalo population of 
our SIDM simulations and the kinematic data of the MW dSphs we 
construct circular velocity curves for the most massive subhaloes 
within 300 kpc halocentric distance for RefP0-3. The dSphs sam- 
ple consists of the 9 galaxies used in Boylan-Kolchin et al. ( 201 lb I: 
Fornax, Leo I, Sculptor, Leo II, Sextans, Carina, Ursa Minor, Canes 
Venatici I and Draco, selected with the criterion Ly > 1O 5 M0. 
The Sagittarius dwarf was removed from this sample since it is 
in the process of interacting strongly with the galactic disc. This 
sample of bright dSphs (plus Sagittarius) is complete within the 
virial radius of the MW (excluding the possibility of undiscovered 
systems hidden in the galactic plane). The Large and Small Mag- 
ellanic Clouds are considerably brighter than the dSphs, and also 
more massive; Boyla n-Kolchin et al.| ^201 lb| > puts a conservative 
lower limit to the maximum rotational velocity of the subhaloes 
that could host these systems: V max > 40 kms -1 . Although in this 
work we do not attempt to find Magellanic Cloud analogues in our 
simulations, we remark that some of the most massive subhaloes 
would actually correspond to these systems rather than to the MW 
dSphs. 

Since the MW dSphs are among the most DM-dominated 
galaxies, having very large dynamical mass-to-light ratios (e.g. 
|Gilmore et a l. 2007), their stars are reliable dynamical tracers of 
their underlying DM subhaloes. Since the available data for dSphs 
provides information only in three dimensions, whereas the phase- 
space distribution function is 6-dimensional, the derived radial 
mass profiles from this data usually carries a degeneracy with the 
anisotropy profile which prevents model-independent constrains on 
the underlying mass profile using a standard Jeans analysis (mak- 
ing for example a cored and a NFW profile equally acceptable 



by the data). Analyses that make use of more information in the 
data can extract more information about the mass profile. This has 
been accomplished not just with consideration of multiple stellar 
components, but also through Schwarzschild modelling {Jardel &] 
Gebhardt 2011 ) and indirectly based on arguments about survival 
of cold substructure within Ursa Minor ( [Kleyna et a l. 2003 1 and 
the distribution of globular clusters in Fornax l Goerdt et al.|2006| >. 
Models that successfully reproduce the observed data tend to have 
roughly the same value of the enclosed mass within the half-light 
radius (e.g. Strigari et al.|2007||Penarrubia et al.|20 08 ; Walke ret al.| 
[20091 |Woif et al.||2010) . In particular, |Wolf et al.|(2010> found 
that the value of the mass enclosed within the 3-dimensional de- 
projected half-light radius rw 2 is accurately given by Mi/a ~ 
3G _1 (of os ) n/a (or equivalently V£ IC (n/ 2 ) « 3 (of os )), where 
(°los) ' s tne luminosity-weighted square of the line-of-sight ve- 
locity dispersion. This approximation is valid as long as the ve- 
locity dispersion profile ai os (R) remains relatively flat out to the 
2-dimensional projected half-light radius. The quality of the kine- 
matic data for the sample of dSphs we have mentioned, guarantees 
the estimation of Mi/ 2 with high accuracy. 

In Figure [8] we show the circular velocity profiles at 2 = 
of the 15 most massive subhaloes (largest peak circular velocity) in 
our simulations for RefP0-3 from top left to bottom right. The black 
symbols with error bars represent the estimates of the circular ve- 
locity at the half-light radius for the 9 MW dSphs of the sample de- 
scribed above (the data was taken from Table I of Wolf £t£]20T0|. 
For the SIDM models (RefPl -3) we show our highest resolution re- 
sults (level 3), but for RefPO the circular velocity curves are taken 
from the level 2 Aq-A-2 simulation (2.8 x 65.8 ~ 184 pc spatial 
resolution) since the level 3 simulation is not converged in that case 
due to the cuspyness of the subhalo density profiles. Convergence 
is easier to achieve for the SIDM models, because they do not de- 
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Figure 8. Circular velocity profiles at z = for the top 15 most massive subhaloes (largest peak circular velocity) of the Aquarius-A halo for the different 
SIDM reference models as given in the legends. The upper left panel shows the standard CDM case, while the bottom panels show two examples of the 
vdSIDM models described in section 



2.1 



Observational estimates of V c i TC (»"i/2 ) for the MW dSphs are shown with black circles with error bars i Walker et al. 
2009 . Wolf et al. 2010 1. All SIDM results are shown at level 3 resolution which is sufficient for convergence due to the subhalo density cores that form in these 
models (see Figures^]and[9j. RefPO is shown at level 2 resolution (2.8 X 65.8 ~ 184 pc spatial resolution), because the CDM subhaloes form cuspy profiles 
which require higher numerical resolution for convergence (see Figure|9}. Clearly, the most massive subhaloes in the CDM model are dynamically inconsistent 
with the MW dSphs, whereas the SIDM subhaloes are consistent with the data. We note that the constant cross section RefPl case is ruled out by different 
observations at the scale of galaxy clusters and is shown here only as a reference. One of the shown subhaloes of RefPl entered already the core-collapse 
regime clearly visible from the circular velocity profiles (see also Figure^Jfor the corresponding steep density profiles). 



velop cuspy profiles, but have constant density cores as shown in 
Figure [7] This convergence is explicitly demonstrated in Figure [9] 
(top panels) where we show the circular velocity curves of the 15 
most massive subhaloes for RefPO (left panel) and RefP3 (right 
panel) at two levels of resolution: level 4 (dashed lines) and level 3 
(solid lines). Clearly vdSIDM subhaloes have essentially converged 
circular velocity profiles, whereas CDM subhaloes are still moving 



towards a more concentrated mass distribution with increasing res- 
olutiorj^] The bottom panels of Figure[9]show the density profiles of 
the five most massive subhaloes at all three resolutions (level 5 as 



Although we do not show the RefPl and RefP2 cases in Figure [5] they 
also show good convergence as the RefP3 case. 
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Figure 9. Top panels: Circular velocity profiles for the 15 most massive subhaloes (largest peak circular velocity) of the Aquarius-A halo (left panel: RefPO, 
right panel: RefP3) at z = for two resolution levels: level 4 (dashed lines) and level 3 (solid lines). Note that the CDM subhalo profiles are not converged 
yet, whereas those of the vdSIDM simulation show reasonable convergence due to the central core that builds up in these cases. Although not shown here, this 
is also true for the RefPl and RefP2 cases. Bottom panels: Density profiles for the 5 most massive subhaloes (left panel: RefPO, right panel: RefP3) at z = 
for three resolution levels: level 5 (dotted), level 4 (dashed lines) and level 3 (solid lines). CDM profiles converge towards cuspy profiles, which are not fully 
resolved with level 3 resolution. The cored profiles of the vdSIDM simulation shows reasonable convergence at level 3. 



dotted lines). The cores are clearly visible in the vdSIDM simula- 
tion along with the better convergence. We therefore conclude that 
the velocity profiles of the SIDM simulations shown in Figure [8] 
are not significantly affected by numerical effects. As for the RefPl 
curves in Figure|7]we note one of the shown subhaloes has already 
entered the core-collapse regime by z=0, which is also clearly visi- 
ble from the circular velocity curves. 

As is evident from Figure [8] the vdSIDM models, which do 
not violate any of the astrophysical constraints, create a subhalo 
population whose most massive systems are dynamically consistent 



with the bright dSps. Although the subhalo cores are not as large 
as in the case of a constant self-scattering cross section (RefPl, top 
right), they are substantial enough to alleviate the tension with the 
data present in the CDM case (RefPO). We emphasise that although 
the most massive subhaloes have essentially the same maximum 
rotational velocities, i.e. the same mass, in the vdSIDM cases as in 
CDM, the vdSIDM subhaloes are much less concentrated and are 
therefore consistent with the observational data. 

Looking at Figure [8] we see that at least three of the dSphs 
(Fornax, Sextans and Canes Venatici I, that have a value of r\i% ~ 
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1 kpc) are inconsistent with the circular velocity profile of the 15 
most massive subhaloes in all our SIDM simulations, particularly 
for the cases with a velocity-dependent cross section. This however, 
does not mean that there are no subhaloes in these simulations that 
are a good match to these galaxies. These subhaloes exist but are 
considerably less massive than the ones we show in Figure [8] and 
are more affected by the limited resolution of our simulations. We 
note that even in the CDM case there seems to be enough subhaloes 
that are consistent with these three galaxies, although a possible 
concern is the high number of massive subhaloes that can be po- 
tentially compared with the dSphs. This is already evident in the 
CDM case: looking at the top leftmost panel of Figure 3 o f|BoylarT| 
|Kolchin et al.| ( |2011b) (corresponding to the Aq-A halo), we see 
that in order to find subhaloes that are consistent with all the 9 MW 
dSphs, it is necessary to consider approximately 30 subhaloes. This 
issue is to be expected for the vdSIDM cases we have explored 
here as well. We believe however that this is a problem that is more 
likely related to the stochastic nature of the abundance of the sub- 
halo population in the Aquarius haloes. For instance, looking at the 
second panel from the top left in Figure 3 of Boylan-Kolchin et al. 
( |2011b) , it seems that Aq-B (another of the MW-like haloes simu- 
lated in the Aquarius project) has approximately 15 subhaloes that 
can be compared with the dSphs, which would clearly eliminate 
the aforementioned problem. Of course, a more detailed analysis 
needs to be done to firmly conclude this. Such analysis would in- 
clude, in addition to possibly re-simulating other Aquarius haloes, 
making a one-to-one matching of the simulated subhalo popula- 
tion with the observed dSphs in a similar way to what was done in 
Boyla n-Kolchin et al.|201 lb] In that work, those subhaloes that are 
considered Magellanic Cloud analogs are removed from the sam- 
ple, and also the subhalo ranking is done according to their mass at 
infall (just before they merged with the host halo), which is a better 
proxy of the potential well that shaped the luminous galaxy than 
the subhalo mass at z = 0. We emphasise that the objective of this 
paper is to present an initial exploration of the vdSIDM models dis- 
cussed in Section [2~l"1 and show explicitly that contrary to the CDM 
case, these vdSIDM models do not predict a number of subhaloes 
which are too concentrated to host the bright dSphs. 



5 SUMMARY AND CONCLUSIONS 

The observed abundance and properties of dwarf galaxies have 
been an enduring challenge for the remarkably successful CDM 
paradigm for more than a decade now. Recent observations with 
high quality kinematic data seem to confirm that at least some of the 
DM-dominated MW dSphs have central density cores instead of the 
high density cusps predicted by CDM (Walker & Penarrubia 201 1 1. 
At the same time, the actual number of observed dwarf galaxies in 
the field (as inferred from the HI ALFALFA survey) seems to be 
in tension with the large abundance of dwarf haloes predicted by 
the CDM model ( |Papastergis et al.|201 l[|Ferrero et al.|201 1| >. This 
seems to be a more serious challenge than the well known excess of 
DM subhaloes compared to the number of observed MW satellites, 
pointed out over a decade ago (Klypin et al. 1999, Moor e~et al.| 
|1999[ >, and that has a viable solution based on the quenching of star 
formation (e.g. |Koposov et al.|2009| l. Another potential "failure" of 
the CDM model has been raised very recently by Boylan- Kolchin| 
|et al.|p01 la[ >, who showed that the subhaloes kinematically associ- 
ated to the MW dSphs seem to be much less concentrated than the 
most massive subhaloes of current CDM simulations. A solution to 
all these challenges based exclusively on processes related to the 



formation an evolution of the luminous galaxies within the dwarf 
DM haloes can not completely be ruled out currently. For instance, 
there are suggestions to create cores in CDM haloes as a response 
of the cuspy profile to energetic feedback. Recent hydrodynamical 
simulations seem to be successful in producing cores if strong feed- 
back mechanisms are invoked (e.g. |Maccio et al.{ 2012, Pon tzen &] 
|Governato|20 11) . It is however important to consider viable alter- 
natives to the CDM model that preserve its success on large scales 
and that, at the same time, are able to solve the small-scale prob- 
lems. 

One of these alternatives is to consider the possibility that DM 
is collisional (SIDM first proposed by Spergel & Steinhardt 2000), 
with self-interactions that are able to isotropise the orbits of DM 
particles and create a core in the inner regions of the haloes. In 
order for this alternative to be feasible, the DM-DM scattering cross 
section should avoid current astrophysical constraints, particularly 
stringent in clusters of galaxies that are consistent with the CDM 
model (for a summary of constraints see Buckley & Fox 2010), and 
at the same time be large enough to produce large density cores in 
dwarf galaxies. 

In this paper we make a first assessment of the viability of 
a class of velocity-dependent SIDM (vdSIDM) models ( |Loeb &| 
|Weiner|2011) using high resolution numerical simulations. More 
specifically we explore a simplified particle physics model where 
the self-scattering between DM particles is set by an attractive 
Yukawa potential mediated by a new gauge boson (either scalar 
or vector) in the dark sector jFeng et al. 2010b, Finkbeine r et al.| 
|201 l||L"oetT & Weiner 201 1 1. This model, naturally predicts a trans- 
fer cross section that depends on velocity in such a way that can 
produce a significant effect in dwarf galaxies and avoid astrophys- 
ical constraints on larger scales (see Section |2~T) . 

We implement an algorithm, based on a Monte Carlo ap- 
proach, to account for DM self-interactions within the framework 
of this SIDM model in the GADGET-3 code for cosmological sim- 
ulations; the description and testing of this algorithm can be found 
in Section [2~2l We use this modified code to run simulations of a 
MW-like halo using the initial conditions of one of the haloes of the 
Aquarius project (Aq-A, Springel et al. 2008 1 exploring two refer- 
ence points within the allowed parameter space of the SIDM model 
for a elastic scattering case (RefP2 and RefP3; see Table[T]and sec- 
tion |3.1[ l. As a reference, we also explored a constant cross section 
SIDM model (RefPl) which is ruled out by observations in galaxy 
clusters (see Sections |2. 1 1 and [3~T| (. We find that for both of these 
cases, the density profile of the main halo remains the same as in 
the standard CDM case outside ~ 1 kpc, whereas inside this ra- 
dius a very small core develops (see Section [3~2| and Figure|4]l. The 
subhalo abundance and radial number density distribution are not 
affected in the vdSIDM models we explored (Figure|6|, which im- 
plies that in the elastic case, these models have the same problem 
as the CDM model regarding the excess of low-mass subhaloes. 
On the other hand, the internal density profile (within ~ 1 kpc) of 
the most massive subhaloes change significantly due to the forma- 
tion of a constant density core in their centre (Figure [7). We show 
that the circular velocity profiles of these subhaloes are consistent 
with the kinematics of the brightest MW dSphs (Figure[8|. Specif- 
ically, and contrary to the CDM case, there are no subhaloes that 
are more concentrated than what is inferred from the kinematics of 
these galaxies. We therefore conclude that the SIDM models ex- 
plored here significantly reduce the tension encountered between 
the observational data and the standard CDM predictions. 

The highest numerical resolution of the simulations presented 
in this work (m p = 4.911 x 10 4 M , e = 120.5 pc) is suffi- 
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cient to obtain converged density and velocity profiles for the more 
massive subhaloes (Figure [9}, but is not enough to reliable trust 
those for smaller mass subhaloes. In a future work, we plan to 
make a more detailed analysis of the subhalo population extending 
to these small-mass systems and use an abundance matching tech- 
nique (based on the mass of the subhaloes at infall, which is better 
correlated with the luminosity of the galaxy) to make a one-to-one 
matching between the simulated subhaloes and the MW dSphs. We 
expect that this will further reduce the tension between the dark 
subhalo population and the observed dSphs, and strengthen the case 
for vdSIDM models. An additional improvement of our algorithm 
is to explore the possibility of inelastic scattering. In particular, in 
models where excited, nearly-degenerate, DM states are possible, 
exothermic collisions could "evaporate" subhaloes by ejecting par- 
ticles during down-scatterings. This could potentially alleviate the 
problem on the overabundance of low mass subhaloes as well. 
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